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ABSTRACT 

The effects of supercooled large droplets (SLD) in 
icing have been an area of much interest in recent years. 
As part of this effort, the assumptions used for ice 
accretion software have been reviewed. A literature 
search was performed to determine advances from other 
areas of research that could be readily incorporated. 
Experimental data in the SLD regime was also 
analyzed. A semi-empirical computational model is 
presented which incorporates first order physical effects 
of large droplet phenomena into icing software. This 
model has been added to the LEWICE software. 
Comparisons are then made to SLD experimental data 
that has been collected to date. Results will be 
presented for the comparison of water collection 
efficiency, ice shape and ice mass. 

NOMENCLATURE 

a p particle acceleration (m/s 2 ) 

A cross-sectional area of particle (nT) 

D drag force (kg*m/s 2 ) 

D a distance between particles (m) 

D k Airfoil leading edge diameter (m) 

d particle diameter (m) 

f drop frequency (Hz) 

g gravitational constant = 9.8 m/s 2 

h film thickness (m) 

L lifting force (kg*m/s“) 

m particle mass (kg) 

s surface wrap distance (m) 

t Time (s) 

V Velocity (m/s) 

V t x-component of velocity (m/s) 

V y-component of velocity (m/s) 

Vol volume (m 3 ) 

x horizontal direction (m) 

x p x-location of water particle (m) 

x first derivative of x-location of water 

particle with respect to time 
(x-component of particle velocity) (m/s) 


x second derivative of x-location of 

water particle with respect to time 
(x-component of particle 
acceleration) (m/s 2 ) 
y vertical direction (m) 

y 0 y-value of the starting locations of 

collection efficiency trajectories (m) 
y p y-location of water particle (m) 

y first derivative of y-location of water 
particle with respect to time 
(y-component of particle velocity) (m/s) 
y second derivative of y-location of 

water particle with respect to time 
(y-component of particle 
acceleration) (m/s ) 

Dimensionless Numbers 

c , coefficient of lift = — 

A p P a V 2 / 2 

c d coefficient of drag = — 

A p p a V 2 ! 2 

,2 

® pP\\P 

Eo Eotvos number = 

a 

f* dimensionless drop frequency = f 

V 

Fr Froude number = 

yjdg 

K Mundo splashing parameter 

5 / 

= <9/ ; Re /4 

K mt Marengo and Tropea parameter 



10000 


La Laplace number = — 
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Oh 

u 

Ohnesorge number = . 

ppad 

Re 

Reynolds number - 

A 

We 

Weber number = ^ — 
a 

Greek Letters 


and impact physics. It will include a review of related 
physics that have been studied in other research areas. 
The second section will describe the modified equations 
including analysis and observations from tests 
performed in the NASA icing research tunnel (IRT). 
The third section will provide comparisons of collection 
efficiency, ice shape and ice mass from those tests with 
the current model. 


LEW1CE 


a angle of attack (degrees) 

P collection efficiency (dimensionless) 

y angle difference between particle 

velocity vector and air flow 
velocity vector (radians) 

8 film thickness (dimensionless = h/d) 

p viscosity (kg/ms) 

V kinematic viscosity of air (m /s) 

p density (kg/m 5 ) 

<7 surface tension (kg/s 2 ) 

0 impact angle (degrees) 

61 empirical splash parameter 

Subscripts 

a air 

1 ice 

n normal direction 

o incoming drop value 

p particle 

s splashed drop value 

t tangential direction 

w water 

x x-dependent 

y y-dependent 

o free-stream property 

INTRODUCTION 

In 1996, the authors presented a report 1 that 
assessed the capabilities of ice accretion software in the 
supercooled large droplet (SLD) regime. This effort 
was performed as a response to the 1994 crash of an 
ATR-72 in Roselawn, IN. 2 It has been speculated that 
accident occurred due to the accumulation of SLD ice 
aft of the deicing boots. Since then, several 
experimental efforts have been made to document SLD 
ice shapes and to investigate the underlying physics. 
This report will present modifications made to the ice 
accretion software, LEWICE, 1 based upon observations 
from those tests. 

The report is divided into three sections. The first 
section will provide a description of the LEWICE ice 
accretion model, with emphasis on particle trajectory 


The computer program, LEWICE, embodies an 
analytical ice accretion model that evaluates the 
thermodynamics of the freezing process that occur 
when supercooled droplets impinge on a body. The 
atmospheric parameters of temperature, pressure, and 
velocity, and the meteorological parameters of liquid 
water content (LWC), droplet diameter, and relative 
humidity are specified and used to determine the shape 
of the ice accretion. The surface of the clean (un-iced) 
geometry is defined by segments joining a set of 
discrete body coordinates. The software consists of four 
major modules. They are 1) the flow field calculation, 

2) the particle trajectory and impingement calculation, 

3) the thermodynamic and ice growth calculation, and 

4) the modification of the current geometry by addition 
of the ice growth. 

LEWICE applies a time-stepping procedure to 
"grow" the ice accretion. Initially, the flow field and 
droplet impingement characteristics are determined for 
the clean geometry. The ice growth rate on each 
segment defining the surface is then determined by 
applying the thermodynamic model. When a time 
increment is specified, this growth rate can be 
transformed into an ice thickness and the body 
coordinates are adjusted to account for the accreted ice. 
This procedure is repeated, beginning with the 
calculation of the flow field about the iced geometry, 
then continued until the desired icing time has been 
reached. 

Collection efficiencies are calculated in LEWICE 
through the use of a particle trajectory analysis. 
Droplets are released from a point in the freestream 
flow and tracked through the flow field using the 
following equations: 

mx p = -D cos y— L sin y+ mg sin a 

my p =-Ds'my+L cosy- mg cos a 

where 


y P -Vy 
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L =C[ 


per 


2 


A 


P 


V = J(x p -V x ) 2 +(y p -V y ) 2 


The initial release point in the x-direction is 
determined by finding an x-location where all of the 
velocities in a vertical sweep are within 0.1 percent of 
the freestream value. The initial release point in the 
y-direction is determined from the angle of attack. The 
initial velocity is assumed to be the terminal droplet 
velocity, which is given by 


c d Ref 


4 gd (Pw-Pq) 
3 V 2 a p a 


where 


Re, 


V,d 


Each droplet is then tracked until it either hits the 
airfoil or reaches the trailing edge. After the first 
trajectory ends, the next particle is released from a 
higher or lower starting point in an attempt to hit the 
surface. This process is continued until there exists at 
least one drop that passes above the airfoil and one that 
passes below. 

Impingement limits are then found using a standard 
bisection search algorithm. Using the coarse limits 
found in the prior step, LEWICE starts a drop halfway 
between these limits. Based upon the end result of that 
trajectory, the next drop is released halfway between 
the current starting point and the starting point of either 
the upper or lower coarse limit. The coarse limit is then 
refined based upon the trajectory results. This process is 
repeated until the starting point of a drop that hits the 
airfoil and the starting location of a drop that misses the 
airfoil is within 1 0 . The bisection is then repeated to 
find the second impingement limit. 

Collection efficiency is then determined by sending 
a user-determined number of trajectories uniformly 
spaced between the impingement limit starting 
locations. The starting and ending locations of these 
trajectories is stored. Collection efficiency is then 
calculated by the following definition: 


This section will present various assumptions used 
in the collection efficiency calculation in LEWICE. 
While other surface physics such as evaporation, rivulet 
flow and film dynamics are important to the water 
impact, much of this review centers on trajectory and 
splashing phenomena. The review will focus on effects 
for larger droplets. A large drop in this context applies 
to any drop size larger than 40 pm, the current upper 
limit in the FAA certification envelope. LEWICE uses 
the following assumptions in the trajectory equations: 


solid particles 

spherical particles 

drops do not break up 

particles do not rotate 

particles have no lift 

particles have no moment 

drag for a stationary sphere applies 

no transient effects of drag 

evaporation of the drop is negligible 

turbulence effects are neglected 

flow is incompressible 

gravity is considered 

drops do not interact with each other 

continuum flow around drop 

all drops that strike the airfoil impinge 


Droplet motion and impact has a wide variety of 
applications and has been studied extensively in several 
research areas. Much of the early work was 
summarized by Clift, Grace and Weber 4 and also by 
Tavlarides et al . ' A summary of more recent research 
was performed by Michaelides. 6 


Third Order Effects 


Third order effects in this context refers to physical 
effects that are negligibly small. Therefore this section 
will review physical effects that can safely be ignored 
in the SLD regime. First, the Basset force models the 
effect of the change in drag as a function of time. This 
term is also called the drag history term. It is used to 
model the flow of gas bubbles, particularly when they 
are immersed in air. This effect is unimportant for cases 
where the viscosity of the surrounding fluid is less than 
the particle viscosity. Since the viscosity of air is much 
less than that of water, this term is safely ignored. 

A lift force for small particles was described by 
Saffman. 7 The Saffman force is often described as a 
"champagne bubble" effect and is more important for 
gas bubbles than for liquid particles. It is also more 
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important for very small particles (smaller than the 
15 pm Appendix C limit) than for the larger particles 
that are being studied. Saffman and Basset forces have 
been calculated in LEWICE and both were found to be 
extremely negligible. 

Turbulence effects have been studied by several 
authors. The works of Mashayek, 8 Marchioli and 
Soldati, 9 and Shin et al. 10 perform DNS calculations of 
turbulence effects on particles, including effects near a 
solid boundary. All of these studies are limited to 
particle sizes much smaller than 15 pm. Particles in the 
Appendix C and SLD matrices are less susceptible to 
turbulence effects. The small effect of turbulence in 
SLD can be seen in the IRT and other icing tunnels. 
Tunnel uniformity depends partly on the droplets being 
dispersed due to turbulence. Some of the difficulty in 
achieving tunnel uniformity can be attribute to this 
factor. Turbulent effects effects have been studiedby 
Bhargava, et al. 11 in their simulation of water droplet 
clouds in the IRT. Compressibility effects are 
considered unimportant unless higher speeds, such as 
those that occur with rotor blades or propellers, are 
encountered. 

Gravitational Effects 

LEWICE models gravity both as a force during the 
trajectory calculation and by using the drop's terminal 
velocity as an initial condition for the drop. As drop 
size increases, the gravitational force will become more 
significant both in terms of this initial velocity vector 
and in the subsequent trajectory path. However, even at 
extremely large drop sizes, gravity has a third order 
effect on collection efficiency. 

This effect was illustrated in the following test 
case. An icing simulation was performed for an MS317 
airfoil with LEWICE using a 10-bin, 236 pm MVD. 
This drop size distribution has the largest MVD of the 
impingement database. 12 The case was computed 
including gravitational force and then repeated with that 
force turned off in LEWICE. This simulates the 
difference in collection efficiency that would occur 
between a wing that is mounted horizontally in a tunnel 
with one that is mounted vertically. As seen in this 
figure, gravity has a negligible effect on water 
collection. The only gravitational differences in 
LEWICE are at the impingement limits. However, 
comparisons with experiments have shown that 
LEWICE overestimates water collection at the 
impingement limits, thus the effect of gravity may be 
even less than shown here. 


MS317 airfoil AOA = 0° 



Figure 1: Effect of Gravity on Collection Efficiency 

Another way to look at the effect of gravity is the 
Froude number. The Froude number is the ratio of 
inertia force to gravitational force and is given by: 



The following plot shows the range of Froude 
numbers that can occur in icing. This plot shows that 
even at extremely high drop sizes, the Froude number is 
close to 1000. A Froude number of 1000 means that 
inertia forces are still 1000 time more important than 
gravitational forces. 



Figure 2: Froude Number Range in Icing Analysis 


Second Order Effects 

Second order effects in this context refers to 
physical effects that can have a minor affect on water 
collection under certain conditions. These effects are 
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modeled by LEW1CE but have been found to be small 
when compared to the inertia force and to the surface 
dynamics of splashing. 

Solid and Spherical Particles 

This assumption states that the drop remains 
spherical and does not deform due to the drop velocity. 
It also assumes that there is no internal flow within the 
drop. The literature reviews cited earlier demonstrate 
that the preferred method of accounting for a change in 
drop shape or internal flow in a drop is to change the 
drag and lift using empirical equations. Beard and 
Pruppacher 13 measured large raindrops falling in air at 
terminal velocity due to gravity and showed that the 
non-spherical shape can be accounted for by using an 
alternate drag model, which at most is 15 percent higher 
than the drag on a sphere. This drag model is shown in 
Figure 3. In addition, their tests showed negligible 
effect on lift and moment. 



Re 

Figure 3: Drag Increase Due to Non-spherical 
Particles 


A test case was then constructed using an MS-317 
airfoil. The following conditions were used: 


a= 2 ° 

LWC= 0.34 g/m 3 
V= 195 mph 
T= -10°C 
MVD=1000 pm 


The first case uses the standard drag model used in 
LEWICE. This drag model uses a correlation of the 
experimental drag values in Figure 1. The second case 
then increased this value 15 percent for the entire 
Reynolds number range. This is a more severe 
requirement than the Beard and Pruppacher model 


where the drag at low Reynolds numbers conforms to 
the standard model. In addition, this case only accounts 
for the drag increase due to deformation and does not 
consider breakup. The actual drag increase will be less 
due to the drop breakup. The collection efficiencies for 
these two cases are shown in Figure 4. As can be seen 
from this figure, the effect of drop deformation on 
collection efficiency is negligible. 



Figure 4: Drag Effect of Non-spherical Drop 


Drop Interaction 

LEWICE assumes that the particle stream is sparse 
enough that interaction effects (such as droplet 
collisions) do not exist or are rare enough to be 
neglected. This assumption allows the solution of each 
drop size in the distribution separately. This assumption 
also neglects the effect of particle density on the airflow 
solution, allowing the airflow to be determined a priori 
(uncoupled). Mulholland et al. 14 studied the effect of 
droplet spacing and showed that for droplet Reynolds 
numbers less than 250, local particle densities greater 
than 0.3 g/m 3 could cause a decrease in drag coefficient. 
However, even for particle densities of 3 g/m 3 , the 
effect on droplet drag was less than 5 percent. 

A second factor that has not been studied 
extensively concerns droplet interactions near the wall 
during a splash event. Smyrnaios et al. 15 studied the 
effect of heavy rain on airfoils and predicted premature 
boundary layer separation for an LWC of 
approximately 3 g/m 3 using a Direct Numerical 
Simulation (DNS) model. A splashing event can result 
in a high density of particles near the wall that may 
interact with the incoming droplets. The effects on the 
boundary layer and the trajectories of incoming 
particles may be taken into account by methods that use 
particle density (by means of droplet frequency) in their 
splashing models. This effect will be discussed in the 
review of splashing effects. 

Drop Breakup 

If a large drop moves at a high enough velocity, it 
can breakup due to shear. Breakup occurs when the 
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drop passes a critical Weber number. Values for this only considers breakup due to the Weber number 

critical Weber number vary widely in the literature. criteria. 

The Weber number is given by: 


We p 


PaV 2 d 

a 


For water drops falling at their terminal velocity, 
the critical Weber number (based on air density) is 
approximately 10. For water drops accelerated by a 
shock wave, a value of 6.5 is given. Krzeczkowski 16 and 
Hsiang 17 each measured droplet breakup for shear 
induced flows and reported values ranging from 10 to 
20. Ibrahim et al. ls provided a more detailed analysis of 
the droplet deformation and breakup using a Taylor 
analogy model. The Weber number of each trajectory 
was output from LEWICE for the case described above 
to investigate this effect. A contour plot of Weber 
number is shown in Figure 5 and shows that the Weber 
number clearly indicates that drop breakup occurs for 
this drop size. 


Corftour 
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Figure 5: Weber Number on 1000 Micron Drop 


Drop breakup is also attributed to the Eotvos 
number, which is given by 


Eo p 



where 


C0nt*ur 

Leoeis 
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4 

& 
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Figure 6: Eotvos number on 1000 micron drop 


The case for a 1000 pm drop clearly shows that 
according to the Weber number criteria, drops will 
breakup before reaching the airfoil. At this point, it is 
unclear how this breakup affects the collection 
efficiency. The smaller drops produced will tend to be 
deflected more, however by the time they reach critical 
Weber number values, they are only 0.1 chord from the 
leading edge even in this extreme example. As most of 
the particle deflection occurs within this region and 
since drops tend to break up into much smaller drops, it 
seems feasible that there is some mass loss that can be 
attributed to this factor. 

An empirical relationship found in Hsiang and 
Faeth 17 was added to LEWICE in order to estimate the 
reduction of collection efficiency due to breakup. In 
their model, droplets will start to break up when the 
critical Weber number is greater than 13. This Weber 
number is based on the air density as defined earlier in 
this report. Since droplet breakup occurs rapidly 
compared to the trajectory time step, breakup is 
considered to be instantaneous. Dai and Faeth 19 
produced some excellent photographs of the breakup 
process using pulsed shadowgraphy and holography. 

Secondary particle size is given by the following 
equation: 


d„ = 6.2 


At 

V Pa J 


- 1 / 

Re„ nd n 


a 


P 


1-2 , -2 
\x +y 


where 


is the acceleration of the particle. 

A critical Eotvos value of 16 or higher is cited for 
drop breakup. The corresponding plot of Eotvos 
number is shown in Figure 6. This shows that although 
the Weber number is high enough to cause breakup, 
Eotvos number is not. The remaining analysis there fore 


Mw 

This correlation can be applied to an Eulerian 
system as well as the Lagrangian tracking system used 
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by LEWICE. However, it would be necessary in a 
Eulerian system to solve coupled sets of equations for 
each drop size generated. In the Lagrangian system, the 
smaller drop size is simply tracked from the breakup 
location. An empirical relationship was chosen to assess 
the importance of breakup to the collection efficiency. 
If breakup is not important then there is no need to 
implement the more complicated droplet deformation 
and breakup (DDB) model described by Ibrahim. 18 

The MVD = 236 pm case shown earlier for the 
MS317 airfoil was the first case chosen to illustrate the 
effect of breakup. Figure 7 shows the case with and 
without breakup. An analysis of the droplet trajectories 
shows that breakup does occur in this case. 



s/c 


Figure 7: Effect of Droplet Breakup on 
Collection Efficiency 


To further illustrate the potential effect of breakup, 
a second example was chosen. In this case, a three 
section multi-element airfoil shown in Fig. 8 was 
chosen since breakup will occur only in regions of high 
velocity gradient. The conditions for this case were 
chosen to represent the most extreme environment 
possible. 



Figure 8: Multi-Element Airfoil Used for Drop 
Breakup Analysis 


For this case, the chord length (stowed) was 
3.96 m, the angle of attack was 8° and the MVD was 
560 pm with a 7-bin distribution. The droplet 
distribution was taken from a draft of the Appendix X 
conditions. 2 " In this case, droplet breakup does affect 
the collection efficiency, especially for the slat. This 


result is expected since droplets near the leading edge 
of a slat will exhibit the highest shear velocities. 


Breakup Effect for Freezing Rain Case - Slat 



Breakup Effect for Freezing Rain Case - Main 



s/c 


Figure 10: Breakup Effect on Main Element 


Breakup Effect for Freezing Rain Case - Flap 



s/c 


Figure 11 : Breakup Effect on Flap 


NASA/TM— 2004-2 12916 


7 


Ice Density 

In the past, LEWICE has always been compared to 
experimental data through the comparison of the ice 
shape profile. These comparisons only use the ice shape 
geometry and not the ice mass. Recently, Potapczuk 19 
measured both ice shape and ice mass on a NACA0012 
airfoil in the IRT. It was shown that ice mass decreased 
as a function of drop size which suggests a mass loss 
due to splashing. However, even though mass 
decreased, many ice shape profiles did not. This result 
demonstrated that the ice density was also decreasing 
with increasing drop size. While the ice mass data 
accumulated so far is insufficient for generating an 
accurate correlation, the following expression was used 
for the ice shape comparisons in this report: 

Pi ~ 920 - 2* d 0 

FIRST ORDER EFFECT: DROP SPLASHING 


The literature search performed for this report 
confirms that the primary assumption LEWICE uses 
that is invalidated for SLD is the assumption that all 
drops that strike the surface impinge, thus neglecting 
splashing and/or bouncing of drops. One of the earliest 
detailed experimental studies was performed by Stow 
and Hadfield 22 who reported on the impact of water 
drops on a dry surface. Macklin and Metaxas 23 reported 
a similar study that also used ethanol and glycerol to 
study the effect of different fluid properties. Jayarante 
and Mason 24 looked at bouncing and splashing of 
raindrops impinging at various angles on dry surfaces 
and films. Wright 25 developed a theoretical splash 
model for raindrops for the purpose of modeling soil 
erosion. 

Harlow and Shannon 26 solved the Navier-Stokes 
equations for the impact of a single drop on a dry 
surface or film. A more recent work was performed by 
Yarin and Weiss 27 28 who proposed a splashing model as 
a type of kinematic discontinuity. Other works include 
Rein 29 who provides a review of several papers, 
including phenomena such as bouncing along with 
splashing and coalescence and Chandra and Avedisian 30 
who documented photographically the droplet structure 
during the deformation process. 

Computationally, a detailed physical model of 
droplet splashing would require solving the Navier- 
Stokes equations for each droplet impact using a 
Volume of Fluid (VOF) model such as that described 
by Hirt. 31 Examples of this type of calculation were 
reported by Trapeaga and Szekely 32 as well as Tan and 
Papadakis. 3 Current computational capabilities usually 
limit this approach to single drop calculations. In a 
typical icing encounter, thousands of droplet impacts 
are recorded per second, making this type of analysis 


prohibitively expensive. An empirical or semi-empirical 
approach is therefore necessary. 

A recent experimental study by Mundo, 
Sommerfeld and Tropea 34 ' 6 examined droplet-wall 
collisions and correlates splashing in terms of Reynolds 
number and Ohnesorge number 


Oh = 


We 

Re 


A 



The Reynolds and Ohnesorge numbers are based 
on the liquid properties and the component of the 
impact velocity normal to the surface. Based on the 
results of their experiment, splashing occurs if the 
factor K=Oh*Re 125 is greater than 57.7. A plot of this 
parameter for drop sizes of 20 and 200 microns is 
shown in Figure 12. 


20 micron, lower surface 
• 20 micron, upper surface 

- - 200 micron, lower surface 

— 200 micron, upper surface 



0.1 0.2 0.3 0.4 0.5 


Figure 12: K-factor for droplet splash 


A small amount of droplet splash is seen in Fig. 5 
even for a 20 micron drop, showing that this 
phenomena will occur at much lower drop sizes than 
droplet break up. This figure also shows that droplet 
splashing is a significant factor in the large drop 
regime. 

The Mundo papers also provide a characterization 
of the size, velocity and direction of the splashed 
particles. Their later references provide an empirical 
splashing model that can be used in Lagrangian 
tracking schemes. The empirical formulas calculate 
splashed drop size, splash velocity, splash angle and 
deposited mass as functions of the incoming 
parameters. The Mundo expressions are given below. 


K = 


r p 3 d 3 V 5V4 

/ w_ y n 


>57.7 for splash 


— = 8.72e _0 0281A: ; 0.05 <-^-<1 
d 0 d a 

n, = 1.676 * 10“ 5 K 2 - 539 ;n < 1000 

A 7 s 
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= 1.337 —1.318— + 2.339 

V t , 0 d Q 
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v„ 


■ = -0.249 - 2.959— + 7.794 


J 


- = n< 


\d 0 j 


Some observations can be derived from these 
expressions. First, for K<77, drops bounce (size out = 
size in). At the splashing onset, roughly half will 
bounce, since half the mass is lost and the size of the 
drops has not changed. At K = 77, all of the drops will 
bounce, as the outbound size has not changed and the 
mass loss is 100 percent. According to this model, the 
breakup of a drop into smaller particles doesn't really 
occur until K>77. Splash mass is 0.15 percent at his 
upper K limit and shows asymptotic behavior. Finally, 
these correlations are only valid up to a drop size of 
150 pm and a drop speed of 18 m/s. This range is much 
lower than that needed for icing analysis. 

However, there are additional problems with this 
approach. If we assume that the splashing measured by 
Mundo scales into the icing regime (higher velocity and 
initial drop size), then maximum mass loss (K = 77) 
would be greatest in Appendix C regime, not SLD. 
Additionally, the correlations given above do not make 
physical sense. For example, this report gives 1000 as 
the maximum number of drops that can be generated 
from a single splash. However, the correlation gives 
only nine drops generated at K = 180, the upper limit of 
his data. Mundo also gives droplet velocity correlations 
that do not conserve momentum at the lower range of 
his K-values. 

Based on these discrepancies, additional reports 
were sought to determine correlations more suited for 
use in icing. Several authors provided empirical models 
for the splashing threshold or for splashed drop size. 
However, a complete model including splashed mass 
loss and splashed velocity was sought so that it could be 
implemented into LEWICE. The following section 
provides a summary of the splashing models found to 
date. 


Summary of Splashing Models 

The parameters needed for an empirical splashing 
model are the splashing threshold, the splashed drop 
size (or drop size distribution), the splashed velocity (or 
a distribution of splash velocities), the splash angle (or 
a distribution of splash angles) and the amount of 
splashed mass. The number of splashed particles is 
needed only if all splashed particles are to be tracked. 


The models presented typically provide these variables 
as a ratio to the incoming parameters. Some of the older 
and less developed splashing models found in the 
literature were not included in this summary. 

Many of these models include the effects of droplet 
frequency, f, and the dimensionless film thickness *. As 
stated earlier, splashed droplets are likely to interact 
with incoming drops. Since this physical effect occurs 
in the experiments reviewed below, it was assumed 
unnecessary to otherwise account for droplet-droplet 
interactions. Droplet frequency can be calculated from 
the liquid water content by assuming that particles are 
uniformly distributed in the freestream. The drop 
frequency is defined as the mass flow rate of water 
divided by the particle mass. The volume of air around 
each drop is assumed spherical with a circular cross- 
section. Therefore, 


. LWC*V*-D l 

4 a 

m 


A 




3 LWCVfDg) 2 

2 p w d y d J 


However, liquid water content (LWC) is the mass 
of water per volume air, thus: 


„ ^ ,3 
Pw , " 


LWC = - 


^Z ) 3 
6 a 


■ = Pw 


l AJ 


Eliminating D a yields the following expression for 
drop frequency: 


3 V 

/ = -- 
2 d 


LWC 

Pw 




The correlations implemented into LEWICE 
calculate frequency from the MVD and the overall 
LWC. If necessary, it could also be calculated from the 
individual droplet spectrum. Film thickness was 
estimated from a correlation provided by Feo: 7 



' LWC' 

[d J 1 

\ Pw y 



The correlations provided in this report were 
developed using similar experimental techniques and 
therefore had similar ranges of applicability. The upper 
limit on drop size and velocity were 340 pm and 30 m/s 
respectively. Droplet frequency and film thicknesses 
were in the range expected for icing encounters except 
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at the lower range. Splash data exists for dry surfaces 
(80) and for film thicknesses of 0.3<8<3. The 
applicability of these models to very thin films is 
unknown. Similarly, droplet frequencies are lower in 
the SLD range due to the lower water contents as well 
as the higher drop sizes. However, the correlations are 
well behaved in these limits and tend toward limiting 
values. 

Lee and Rvou Model 

Lee and Ryou 38 proposed the following statistical 
model: 



d 3 
' V ~ 2 


LWC 


P 


w J 


pdf 


r d s ) 
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(d s /d a ) 

b - 1 

exp 


(d s /d 0 ) 
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where 



= -0.993 + \.160 o -1.56 6> 0 2 + 0.496>o 

n,o 

— = 0.2 + 0.97W(0,l) 
m Q 

n s =0.1 87We 0 - 4.45 


£ = 2.71-9.25*10 (fe 


61, =0.21-7.69*10 ~ 5 We 


Pdf 


f V ^ 

r n,s 




r V /V 

r n,s' r n,o 




V ®o ) 


exp 


f V IV 

r n,s' r n,o 


'A 


o v 


where 


b =2.1 for 9, < 50° 


b v = 1.1 + 0.02 6j for 6j > 50° 





where RN(0,1) is a random number between 0 and 1. 
This approach requires a large number of trajectories to 
be calculated in order to achieve a statistical sampling. 
This occurs because the splashed mass and splashed 
drop size are calculated using a random number. The 
underlying assumption is that the phenomena is chaotic. 

Stanton and Rutland Model 


Q v = O.158e°'° 1700 

6> = 65.4 + 0.266 6> 

V, = V„ cot 0 s 

Marengo and Tropea Model 

Marengo and Tropea 41 used a similar form for drop 
size and velocity components: 


Cossali et al. 39 summarized a number of splashing 
models and made comparisons to their data. One of 
these models was derived by Stanton and Rutland 40 
using a complex correlation from the data of Yarin and 
Weiss. 28 They used a distribution of splashed droplets 
rather than a single splashed drop in their model. This 
approach would also require calculation of a large 
number of splashed droplet trajectories since each 
incoming drop would produce a distribution of splashed 
particles. The correlation is given by the following 
equations: 


-^= (0.25 +0.238 S)- (0.22 +1.28 5\k mt -K mt c ) 
d„ 


K mt =K /5 f 


8/ *- 6/5 


HO" 


^ = (0.056 + 0.057£)+ 0.38 {K mt -K mt J 

" t n 


= ( 0.311 - 0 . 077 ^)- ( 0.09 + 0 . 24 5){K mt - K mt c ) 

" n,o 
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(0.3 63 + 0.242 5){K mt - K mlc f- 92S ~ l 52 w 


n s = 641.8 + 


640 . 8 (K mt K mtc ) o (o.675+o.o36<y)(^„ -k„ c ) 


Trujillo Model 


Samenfink Model 

Finally, Samenfink et al. 44 used the following 
expressions: 

— = 1 -0.03454S 0 ' 175 ^' 1239 Ta 0 ' 265 

d n 


Trujillo et al. 42 used Mundo's drop size expression — -0 OS 0 145 -0 ' 3384 # 0 ' 2938 ( y _003113 / jfl 01157 

but used different forms for the other variables. Their V 0 

V O 

correlation can be expressed by: 


V ts 
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1/ * -3 / 


1 / „ * 


- 3 / 


K /2 f ^~K/2f c 


22 


0.0437 


K 


( v ^ 2 

r X,0 


\ V v,oj 


-K r 


-44.92 


Schmehl Model 


Schmehl et al. assumed that the splash velocity 
was 60 percent of the incoming velocity and used the 
following expressions for the remaining terms: 


In — L = — 2 — 
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In this model, the splashing parameter S has the same 
definition as used in the Schmehl model and splashing 
occurs for S > 1 . 

COLLECTION EFFICIENCY RESULTS 

By knowing the splashing parameters, a feature can 
be added to LEWICE to track the trajectories of the 
splashed particles and the trajectories of particles after 
breakup. This process was described in a recent report 
by Rutkowski, et al. 45 Mass loss due to splashing and 
subsequent reimpingement are also calculated. Since it 
is necessary to extrapolate from the experimental data 
for icing encounters, care must be taken, especially in 
high velocity impacts. The Trujillo model and the 
Samenfink model were programmed into LEWICE to 
assess the effects of splashing on water collection and 
ice accretion. These two models were chosen as the 
correlations tended to exhibit asymptotic behavior at 
higher velocities. 

The MS-317 airfoil was chosen to demonstrate the 
effect of these splashing models on water collection. 
The MS-317 airfoil has been used in several IRT entries 
for measuring collection efficiency. 1246 The conditions 
chosen were 0° and 8° angle of attack, and drop sizes of 
21 pm, 92 pm, and 236 pm. The 21 pm cases were 
chosen to show that although the models predict 
splashing in Appendix C conditions, very little mass is 
lost. The tunnel speed for all cases was 175 mph 
(152 kts, 78 m/s). 
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As can be seen in Figures 13 and 14, the collection 
efficiencies are very similar for the splashing and non- 
splashing cases. Slightly more water is lost using the 
Trujillo model than the Samenfink model, but the two 
results are comparable. 


MS-317 airfoil AOA = 0° 



s/c 

Figure 13: Collection Efficiency for 
MVD = 21pm and AOA = 0° 


MS-317 airfoil AOA = 8° 

0.9 
0.8 
0.7 
0.6 



s/c 


Figure 14: Collection Efficiency for 
MVD = 21pm and AOA = 8° 


More mass loss is predicted for the 92pm cases 
shown in Figures 15 and 16. The models predict the 
mass loss at the leading edge very well, but do not 
capture the mass loss at the impingement limits. This 
result is to be expected from the empirical models since 
the splashing parameters K and S are calculated using 
the normal component of the impact velocity. 


MS-317 airfoil AOA = 0° 



Figure 15: Collection Efficiency for MVD = 
92pm and AOA = 0° 


MS-317 airfoil AOA = 8° 



Figure 16: Collection Efficiency for 
MVD = 92pm and AOA = 8° 


Figure 17 shows a plot of the impact velocity and 
the impact velocity components for a 241 pm drop. The 
241 pm drop size was the largest drop used in the 
92 pm distribution. The impact velocity is almost the 
same at each chordwise location. This occurs because 
the drop size is so large its trajectory is almost ballistic. 
The normal component of the velocity is largest at the 
leading edge and decreases chordwise. Splashing will 
therefore have the largest effect at the leading edge 
according to these models. However, the splashing 
parameters in this regime are much greater than the 
upper limits of the experimental data. This suggests that 
while impingement near the leading edge may still be 
modeled using these equations, a new model is needed 
at the impingement limits. 


no splashing 
■ ■ ■ Trujillo model 
— — Samenfink model 
□ Experiment 
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ICE SHAPE AND MASS RESULTS 

The Trujillo model was then used along with the 
ice density correlation to generate ice shapes on a 
NACA0012 airfoil. The dataset chosen was reported 
last year by Potapczuk. 21 This data was chosen as it also 
contained measurements of ice mass as well as ice 
shape. The conditions are shown in the following table. 


-0.1 -0.05 0 0.05 0.1 

s/c 

Figure 17: Impact velocity for a 241pm Drop 


Figures 18 and 19 show the collection efficiency 
results for the 236 pm case. This case represents the 
largest drop size for which collection efficiency data 
exists. The Trujillo model exhibits a similar behavior at 
236 pm as it did for the 92 pm case, demonstrating that 
this model can be extrapolated somewhat from the 
experimental limits. The Samenfink model continues to 
show much higher mass loss at the higher drop size. 
This result is not considered realistic. 

MS-317 airfoil AOA = 0° 


Table 2: Inputs for Ice Shape Cases 


Case 

MVD 

(pm) 

LWC 

(g/m 3 ) 

V (m/s) 

T W (°C) 

time 

(sec) 

Run 

1-17 

40 

1.02 

77 

-19.3 

588 

Run 

1-1 

70 

0.91 

51 

-19.6 

804 

Run 

1-4 

160 

1.5 

52 

-19.5 

300 

Run 

1-22 

40 

1.02 

77 

-19.3 

576 

Run 

1-23 

70 

0.65 

77 

-19.2 

714 

Run 

1-26 

160 

1.04 

77 

-19.2 

336 




Figure 18: Collection Efficiency for 
MVD = 236pm and AOA = 0° 


MS-317 airfoil AOA = 8° 



Figure 19: Collection Efficiency for 
MVD = 236pm and AOA = 8° 


The first case shows the comparison for Run 1-17 
and Run 1-22, the two 40 pm drop size cases. The two 
conditions are very nearly identical and show a similar 
result. The comparison shows a slight decrease in the 
overall shape. Additionally, the shape is closer to a rime 
shape like the experiment. It should be noted that 
LEWICE 3.0 will generate the same ice shape as 
LEWICE 2.0 if the splashing model is deactivated. 



Figure 20: Ice Shape Comparison for Run 1-17 


Run 1-22 



Figure 21 : Ice Shape Comparison for Run 1 -22 
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Run 1-26 


Figures 22 and 23 show the results for the 70 (tm 
cases. The conditions are scaled to produce the same 
accumulation parameter. These cases show very little 
difference between LEWICE 3.0 and the previous 
model. The mass loss due to splashing is nearly equal to 
the reduction in ice density, resulting in a similar shape. 


Run 1-1 



Figure 22: Ice Shape Comparison for Run 1-1 



x/c 

Figure 23: Ice Shape Comparison for Run 1-23 


Figures 24 and 25 show the ice shape comparisons 
for the 160 pm cases. Both cases produce a smaller and 
more streamlined shape than LEWICE 2.0 produced 
without splashing. The effect is more pronounced for 
Run 1-26 because the higher velocity will produce more 
splashing and hence mass loss. 



Figure 25: Ice Shape Comparison for Run 1-26 


Finally, Figure 26 shows the prediction of ice mass 
for these conditions. The mass predicted by LEWICE 
3.0 is closer to the experimental data, but this is because 
the ice density correlation was determined from these 
results. More cases are needed to determine if this result 
is real or a by-product of the correlation. More 
comparisons, including glaze ice shapes, are needed to 
fully assess the splashing model. 
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Figure 26: Ice Mass Comparison for Example Cases 


Run 1-4 



X/C 

Figure 24: Ice Shape Comparison for Run 1-4 


CONCLUSIONS 

A review of the physical effects of large droplet 
phenomena confirmed that droplet splashing is the 
primary physics that is not modeled in icing software. 
Several other phenomena such as drop breakup, drop- 
drop interactions and changes in drag, lift and ice 
density were also covered. Empirical models for these 
phenomena were presented. The resulting effect on 
water collection, ice shape and ice mass for LEWICE 
were presented. 
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The results show that the current empirical 
splashing models could only account for splashing 
effects near the leading edge. The models, all of which 
are based on the velocity component normal to the 
surface, could not explain the measured decrease in 
water collection near the impingement limits. However, 
it was necessary to extrapolate values from the 
empirical data for this effort. It is hypothesized that 
droplets that impact at high total velocities (> 30 m/s) 
but low impact angles (< 20°) may rebound in a manner 
different than the current experimental data at the lower 
velocities. 

Ice shape profiles show similar results to previous 
LEWICE outputs since an ice density correlation was 
added to counteract the loss in mass due to splashing. 
The ice mass results for the current model are closer to 
the experimentally measured values, but the same 
dataset that was used to generate the correlation was 
used to assess the results. 
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